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Abstract. The theory of slow manifolds is an important tool in the study of 
deterministic dynamical systems, giving a practical method by which to reduce 
the number of relevant degrees of freedom in a model, thereby often resulting in 
a considerable simplification. In this article we demonstrate how the same basic 
methodology may also be applied to stochastic dynamical systems, by examining 
the behaviour of trajectories conditioned on the event that they do not depart the 
slow manifold. We apply the method to two models: one from ecology and one 
from epidemiology, achieving a reduction in model dimension and illustrating the high 
quality of the analytical approximations. 
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1. Introduction 

It is possibly only a slight exaggeration to say that of all the mathematical models 
we can dream of, there are only two kinds which are straightforward to solve: those 
which are linear, and those which are one-dimensional. This aphorism holds equally for 
stochastic dynamical systems as it does for their deterministic counterparts. Much of 
applied mathematics and theoretical physics is devoted to the delicate art of taking high- 
dimensional or non-linear problems of interest and finding appropriate approximation 
schema with which to reduce their apparent difficultly. 

The theory of 'slow' or 'centre' manifolds is an example of just such a scheme and 
one which is well developed for deterministic dynamical systems [T]. In many models of 
interest there exists a separation of time scales between some quantities which relax very 
quickly to an essentially static value, and others which change more slowly and can be 
sensitive to perturbations. The term 'slow manifold' describes the space in which these 
slower quantities vary, after any fast initial transient has died out. Restricting attention 
to this space offers an effective method by which to remove the so-called fast degrees of 
freedom, thereby reducing the dimension and simplifying the system. In this article we 
will show how the principles of slow manifold theory from deterministic systems may be 
put to good use in stochastic systems as a tool to separate time scales and reduce the 
dimensionality of the model. 

The goal of removing fast degrees of freedom from stochastic systems has received 
significant attention and a variety of approximation methods have previously been 
proposed pHTT]. though a simple and generally applicable theory is yet to emerge. 
The existing literature on the subject can be coarsely split according to the framework 
within which the stochastic system is represented. Some authors have focused on master 
or Fokker-Planck type equations which govern the time evolution of the probability 
distribution of system states [I2]. Within this formalism, a reduction of dimension 
can be achieved through the application of projection operators, as illustrated by 
Gardiner |4], and developed by others [5l|6]. Unfortunately, the details of the method 
can be somewhat cumbersome, especially in cases where the fast degrees of freedom 
are not parallel to the variables with which the system is described. The alternative 
formalism of stochastic differential equations (SDEs) provides a perhaps more intuitive 
approach. The equations describe the evolution of trajectories themselves and bear 
a useful resemblance to ordinary differential equations, which aids physical reasoning. 
For this reason we choose the SDE formulation as the basis for our work. However, 
it should be pointed out that such systems cannot be correctly interpreted without 
specifying the choice of stochastic calculus, and that certain manipulations are not 
entirely straightforward |13] . 

Perhaps the simplest implementation of time scale separation in the SDE setting is 
the treatment of 'direct adiabatic elimination' presented in [7]. It will be instructive to 
review the salient points. The procedure mirrors closely that of slow manifold theory: 
the variables associated with the fast dynamics are assumed to be stationary, from which 
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a function describing the slow manifold may be determined. If the fast variables are 
subject to noise, then the procedure will yield an expression for the slow manifold which 
contains noise variables; unfortunately this can lead to ill-defined terms [7]. The method 
can, however, be usefully applied in systems where the entire system is linear in the fast 
variables. The more rigorously derived and well-known Haken slaving principle [2l|3l[Tl] , 
and several other related methods [HHTO] are developed along similar lines of reasoning 
and suffer from the same complications which arise from specifying the slow manifold in 
a stochastic sense. A notable exception is [15], which deals with a particular model in 
which there is a true centre manifold (that is, a surface on which there is no deterministic 
flow), and applies a novel method in which stochastic perturbations away from the 
manifold are assumed to instantaneously relax along the deterministic trajectories. We 
have applied a similar procedure in several previous works [IIl[16l[l7] , and it is this idea 
which will form the basis of the method we develop here. 

The core of our approach is to examine the behaviour of a stochastic system in 
the SDE framework under the condition that its trajectories are confined to the slow 
manifold of the deterministic version of the system. Because we use a static description 
of the slow manifold, the method is applicable to a broader range of systems than the 
direct elimination procedure or the Haken slaving principle. Moreover, the procedure is 
mathematically explicit and straightforward to apply. One also gains a sense of physical 
intuition as to the behaviour of the system, which is arguably not present in the master 
equation setting. 

We develop and demonstrate the method with the aid of a pair of motivating 
examples, which we select from the burgeoning field of research into the effects of 
demographic noise in individual-based models (IBMs). Models of this type have recently 
become very popular in several areas of physics and applied mathematics as tools to 
study complex emergent phenomena. In the limit of large system size, the mesoscopic 
behaviour of a given IBM is well described by an SDE system derived from the rules of 
the model (see [T^ and the appendices of this paper for details). It is to these SDKs 
that we will apply our method. 

We will also show how the slow manifold approximation can be effectively combined 
with other stochastic approximation techniques. The linear noise approximation (LNA, 
or van Kampen expansion [18]) has found favour amongst theoreticians studying IBMs. 
In a nutshell, the LNA provides a mesoscopic description of the system in terms of a 
linear SDE. The price paid for this simplification is that the theory only applies in the 
neighbourhood of an attractive fixed point of the noise-free version of the model. In 
the context of the LNA, slow manifolds can be a malign presence. If some eigenvalues 
of the approximate linear SDE are close to zero, then a small stochastic fluctuation in 
the direction of the corresponding eigenvector can carry the system very far from the 
steady state into regions in which the true non-linear nature of the model is important. 
Moreover, a large separation between eigenvalues can in some situations lead to the 
numerical evaluation of theoretical solutions becoming ill-conditioned. Both of these 
effects can lead to a very poor agreement between stochastic simulations and the LNA 
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theory. 

In the next section we develop our method with the aid of a simple illustrative 
example from ecology, before providing a general formulation. The results from this 
example demonstrate the success of the approximation scheme even in regimes where 
the fixed point is weakly unstable; this addresses the first difficulty with slow manifolds 
and the LNA identified above. In section three we go on to apply the general formulation 
of the method to an epidemiological model with seasonal forcing. This model has been 
identified as suffering from the technical numerical difficulties associated with a large 
separation between eigenvalues [19]. We show how our method may be used in tandem 
with the LNA to provide a very good approximation to results coming from stochastic 
simulations. 



2. Method 

2.1. Motivation 

We begin by recapping the basics of slow manifolds in deterministic systems. Consider 
the ordinary differential equation 

(1) 

where x is an n-dimensional vector describing the state of the system, and A is an n- 
dimensional vector-valued function of x. As is well known, the behaviour of the system 
in the neighbourhood of a fixed point x* is described by the linearization of A around 
that point. Define the Jacobian matrix J with entries 

Then for x close to x*, the time evolution of ^ = a; — a;* obeys 

Further insight is gained by considering the eigenvalues and eigenvectors of J. From (j3]), 
we learn that if is a left eigenvector of J with eigenvalue A, then errors in the direction 
of V will grow exponentially if A > and shrink exponentially if A < 0. If, on the other 
hand, A = then we do not know what effect perturbations in the direction of v will 
have on the long-term behaviour of the system. To answer this question would require a 
more detailed non-linear analysis, which is likely to be very difficult in a general system 
with several degrees of freedom. Slow manifold theory offers a way to make progress in 
this case by reducing the dimension of the model. 

The basic observation behind the theory is as follows: since perturbations in 
the direction of stable/unstable eigenvectors will shrink/grow exponentially, the only 
trajectories whose behaviour is in question are those which are tangential to the span of 
the eigenvectors with eigenvalue zero. Very often, this set of eigenvectors has many fewer 
members than there are degrees of freedom in the original system, and thus restricting 
attention to this subspace achieves a considerable reduction in dimensionality. 



Stochastic dynamics on slow manifolds 



5 



Slow manifolds are also of great practical use when no eigenvalues are precisely 
zero, but there is a separation of time-scales. For example, suppose a stable fixed point 
X* has associated eigenvalues satisfying Re[Ai] < . . . < Re[Am] ^ Re[Arrt+i] < . . . < 0. 
Perturbations in the direction of eigenvectors Vi , . . . , will decay extremely rapidly 
in comparison with those in the directions of fm+i , • • • ,Vn- For practical purposes, 
the 'slow' manifold of trajectories tangent to these less stable eigenvectors defines an 
(n — m)-dimensional system which will provide a good qualitative approximation to 
the behaviour of the larger system, as perturbations away from this manifold will very 
quickly collapse. 

Our goal in this article is to apply the basic ideas of slow manifolds to stochastic 
systems. Our starting point will be the SDK 

^ = A{x) + n{t), (4) 

where x and A are as in ([T]), and 77 is a vector of Gaussian white- noise variables with 
zero mean and correlations 

(r/,(tH(t')) =£5(t-t')%(^)- (5) 

Here (■ ■ ■) denotes averaging over the noise, £ is a small parameter governing the strength 
of the noise, 5 is a matrix-valued function of the system state, and the SDE (jlj) is to 
be interpreted in the Ito sense. The exact form of A and B will be determined by the 
application. 

We are interested in the case where the deterministic {e = 0) system exhibits a 
slow manifold. How will the stochastic system behave if we confine its trajectories to 
this slow manifold? We develop the theory with the aid of a specific example. 

2.2. Illustrative example 

To illustrate our method, we explore the behaviour of a simple ecological model of two 
interacting populations, labelled X and Y. Individuals of both populations reproduce 
with rate one, and there is a small probability /i of the offspring mutating from one type 
to the other. The organisms also prey on each other with rate e and a slight preference 
p for prey of the opposite type. The model may be written in the traditional notation 
of chemical reaction systems: 

Reproduction: X X + X , Y Y + Y 

Mutation: X — % X + Y , Y — % X + Y 

£{l/2+p) £{l/2+p) 

Predation: X + Y > X , X + Y > Y 

e{l/2-p) e{l/2-p) 

Cannibalism: X + X > X , Y + Y > Y . (6) 
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Here arrows denote possible reactions and the values above them are the rate constants. 
Writing nx and ny for the number of individuals in each population, the model may 
be mathematically formulated as a master equation describing the time evolution of the 
probability distribution P(nx,^y); see (I A. 2 1) in Appendix A. Stochastic simulations of 
the model can be performed efficiently using the Gillespie algorithm [20] . 

When the predation rate e is small, the population may grow very large. Performing 
an expansion of the master equation in the limit of small e yields an effective description 
of the system in terms of an SDE for the scaled variables x = enx and y = eny- We 
give the details in Appendix A. For the present model, we find the following pair of 
equations: 

dx (\ \ 

— = x-^{x-y)-xl-{x + y)-p{x-y)j+ri^{t), 

(7) 

^ = y + fi{x-y)-y Q(x + y) + p{x - y)^ + riy{t) , 

where r]x and r]y have the correlation structure specified in (j5]), with 

X + \x{x + y) - {px + ii){x - y) \ 

\ y + \y{x + y) + {py + ii){x - y) j' 

We begin by examining the deterministic system found by putting e = 0. There 
is a trivial fixed point at x* = , ?/* = 0, representing the extinct state, which is 
always unstable. There is a second fixed point at = 1 , = 1, representing equal 
coexistence of the two populations. This state is stable when p < ^. If p is raised above 
/X, a supercritical pitchfork bifurcation occurs, with the equal coexistence fixed point 
becoming unstable and giving rise to a symmetric pair of stable fixed points in which 
one species dominates the other. The new fixed points have coordinates 



. _ l-2^± v/(l-/x/p)(l-2^) , _ l-2^TA/(l-y^/p)(l-2^) 
^ ~ l-2p ' ^ ~ l-2p • ^ ^ 

We are interested in examining the effect of noise near this transition. 

The eigenvalues of the Jacobian at the coexistence state are —1 and A = 2{p — /i), 
with corresponding eigenvectors (1,1) and (1,-1). If |A| ^ 1 then we have a slow 
manifold in the direction of {x — y), meaning that perturbations to the balance of 
populations evolve very slowly. Formally, the slow manifold is defined by the collection 
of trajectories which are tangent to the slow eigenvector at the fixed point, although in 
practice there is unlikely to be a closed analytic expression for this surface. To make 
progress we approximate the slow manifold by the surface on which the rate of change 
in the direction of the fast eigenvector is zero (known as the nullcline). In the present 
model, the nullcline of the fast direction (x + y) is given by dx/dt + dy/dt = which 
yields the hyperbola 

{x + y)-hx + yy +p{x-yf = 0. (10) 
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Figure 1. These plots show the behaviour of our example illustrative ecological model 
either side of the pitchfork bifurcation. The fixed points are shown as red circles, 
and the dashed red line is the nuUcline dx/dt + dy/dt = 0, given in (fTO]) . The grey 
arrows show trajectories of the deterministic system, while the black line traces out 
the trajectory of a single (short) stochastic simulation of the individual-based model, 
starting close to the origin. The parameters are e = 0.005 and p = 0.3 in both plots, 
while n = 0.35 on the left and /x = 0.25 on the right. 



The two plots in figure [T] capture the typical behaviour of the model for parameters 
either side of the transition. 

The SDE system ([71) is two-dimensional, non-linear and has noise correlations which 
depend on the state of the system. These factors combine to make the theoretical 
analysis of the model very difficult. The situation is not hopeless, however, as it is 
clearly visible in figure [1] that the system state does not typically stray very far from 
one-dimensional subspace defined by the nuUcline of the fast variable (x + y). We intend 
to exploit this fact to produce an 'effective' one-dimensional description of the model. 
The plan of attack is as follows: first we will make a coordinate transform to separate the 
fast and slow variables; then we will examine the behaviour of the slow variable under 
the assumption that the fast variable relaxes instantaneously to its nuUcline value. 

We introduce w = x + y and z = x — y, so that 






In the new coordinates the nuUcline is described by the equation w — w"^ /2 + pz^ = 0, 
and (I7j) becomes 

dw 1 



(12) 
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To determine the correlation structure of the new noise variables and rj^, we apply a 
general result on Gaussian random variables. Suppose that a vector of Gaussian random 
variables 77 has correlation matrix B, and that <^ = Vt] for some matrix V. What is the 
correlation matrix of Well, 

(00) = ( Yl VrkVkVjiVi) = ^ik{VkVi)ViJ = B,, , (13) 

where B = VBV^ . In the present case, the matrices B and V are given in equations 
(IH]) and ( fTTI) . respectively. We thus find the following correlation matrix for rjuj and r^^: 

^^f w + lw^-pz^ z(l -2/i + w(i -p)) \ 
y z{l - 2ij + - p)) w + ^w^-pz^ J' 

Notice that whilst the original noise variables 7]^ and rjy were independent (the off- 
diagonal entries of B were zero), the noise variables in the new coordinates are correlated 
with each other. 

To enforce the assumed separation of time-scales between w and z, we impose the 
following conditions: 



w = 1 + ^/1 + 2pz^ and r/^(t) = . (15) 

The first of these sets w to its value on the nullcline for a given z, whilst the 
second removes the possibility of any noise-induced fluctuations. What effect do these 
constraints have on the evolution of First, we may substitute w = 1 + a/1 + 2pz'^ 
into f|T2|) to remove the dependence on w, thus 
dz 

-^^=f{z)+v.{t), (16) 

where 

f{z)=z(^l-2ii- Q-p^ (l + ^l + 2pz^)^ . (17) 

Second, we must determine the effect of the conditions on the noise variable 77^. Since 
rjw and rjz are correlated, imposing 1]^, = will alter the statistical distribution of r/^. 

Again we apply a general result on correlated Gaussian random variables (see 
Appendix B of [K]). Suppose that a collection of Gaussian random variables (771, ... , rjn) 
has correlation matrix B. Let B be the correlation matrix of (772, . . . ,7/„) conditioned 
on the event that rji = 0. Then B and B are related by 

lB'% = [B-%, for all 7, j = 2, . . . , 7z . (18) 

This fact straightforward to prove by integration of the Gaussian probability density 
function. In particular, if 77 = 2 then the variance of 772 conditioned on 7/1 = is 

B = B,,-^^. (19) 

Applying formula f[TI?l) to the correlation matrix found in ([T^ . we obtain 

{7],{t)r],{t'))=e6{t-t')g{z), (20) 
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Figure 2. The stationary distribution oi z = x — y Sts predicted by the reduced- 
dimension model (orange curve) and measured from a single long simulation run of the 
individual-based model (black histogram). The parameters are the same as those in 
figure [TJ The theoretical prediction was obtained by numerically integrating (l22t . with 
parameters taken from figure [1] whilst 2000 data points were collected from simulation 
run at time intervals of 100. 



where the noise strength g is given by 

1 




g{z) =[w+^-w'-pz']\l-{z \ ] 1, (21) 



and where of course w = 1 + a/1 + 2pz'^. 

Equations (fT6ll and (l20l) together define a one-dimensional stochastic differential 
equation. Although it may look a little complicated, being one-dimensional means that 
most questions of interest about this system can be answered by standard methods [T2] . 
For example, the long-time average behaviour of the model is captured in the stationary 
distribution of z, which has the following explicit form: 

P(^) = J.^exp f- r ^-^dz'] . (22) 



9[z) \e g{z') 
In figure [2] we compare the analytical prediction of fl22|) with a histogram of the z- 
coordinate of the sample points of stochastic simulations taken from figure [H Clearly, 
the reduced one- dimensional model provides a very good fit to the data coming from 
the individual-based simulation. It should also be pointed out that although we have 
developed the theory based on the local behaviour around the coexistence fixed point 
(1, 1), the approximation remains successful even in the unstable regime. 

2.3. General formulation 

We close this section by providing a description of the method for an arbitrary n- 
dimensional SDE 

d T* 

— = A{x) + 'n{t), (23) 
with noise correlations 

'MtHit'))=e6it-t')\B{x)] , (24) 

/ L Jij 
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where i,j = l,...n. 

Suppose we are interested in behaviour around a fixed point x*. Let J be the 
Jacobian of A at that point and write Ai , . . . , A„ for its eigenvalues. Suppose further 
that Ai is non-degenerate, real, large and negative, thus its associated eigenvector Vi 
represents a very stable direction. We aim to eliminate fiuctuations in this direction to 
produce a reduced-dimension model. To apply our method, we first make a change of 
variables from the n-vector a; to a single variable w and an (n — l)-vector z, via the 
coordinate transformation 

=Vx. (25) 
It is possible to choose V so that 

J = V.W= ( r I , (26) 




where L is an {n—l)x[n — 1) matrix. The first row of V must be vf, and therefore near 
the fixed point x*, the variable w describes the distance from the slow manifold along 
the fast direction vi, while the remaining n — 1 slow degrees of freedom are captured by 
z = (Z2, . . . 



)'^. In the new coordinates the SDE becomes 

d I'm 
dt \z 



"f-^^VAfv-^iZ'W^m, (27) 



where 



Q{t)Q{t')'f=e6{t-t')[B{w,z) 



(28) 



and B = VBV^. We wish to constrain w to its nuUcline, which is defined by the 
equation 

v^a(v~'(^^^ =0. (29) 

We will assume that we may unambiguously describe the nuUcline by a known function 
w = 0{z); in the illustrative example, 6{z) was given by (flSi) . a branch of a hyperbola. 
In general however, it may have a more complicated form. 

To enforce the assumed separation of time-scales between w and the other variables, 
we impose the following conditions: 

w = 9{z) and Ci(^)=0. (30) 

For the remaining variables, we have 

^ = A{z) + at) , (31) 

where 

(O(t)O(t')) =e6{t- t') [Biz)] ^ , z,j=2,...n. (32) 
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The drift vector A{z) and diffusion matrix B{z) are derived from A and B as follows. 
For i, j = 2, . . . ,n 



Aiz) 



VA V 



z 



and 



'b{z)-^_ 









B{9{z) 



(33) 



(34) 



The general solution for the new noise covarience matrix can then be shown to be 

Biz) = B22{z) - B2iiz)Bi,\z)B^2iz), (35) 
where B are matrices of the partitioned B matrix; 



B 



Bn{z) Bu{z) 
B2i(z) B22iz) 



(36) 



with -611(2) a 1 X 1 matrix and -622(2) an [n — 1) x [n — 1) matrix (see Appendix B 
of [IS])- Equations (jHID and (15^ describe a reduced-dimension stochastic system in 
which the fast direction associated with the eigenvalue Ai has been eliminated. 



3. Application: seasonally forced epidemics 
3.1. Model definition and deterministic treatment 

The SEIR model is a simplified epidemiological model describing the spread of a 
disease through a population [21]. Members of the population may be in one of four 
states: susceptible (S), exposed (E), infectious (/) and recovered (R). The susceptible 
individuals come into contact with the infected and become exposed with infection rate 
/3(t), which may vary with time. Those exposed to the disease then become infectious 
with a rate of disease onset a. Finally, the infectious recover with an average rate of 
7. In addition to these disease dynamics, there is a constant birth and death rate /i; 
it is traditional to hold the population size constant by treating death and birth as a 
single process whereby an individual returns to the susceptible state. As in the earlier 
illustrative model, the dynamics may be conveniently summarized using the notation of 
chemical reactions: 

Infection: S + I E + I Incubation: E I (37) 

Recovery: I ^ R Death/Birth: E,I,R-^ S. 

We write ns ,nE ,ni ,nR for the number of individuals in states S, E, I and R, 
respectively. The total population size is then given hj N = ns + nE + nj + n^., which 
does not vary, meaning that there are three degrees of freedom. With just a slight abuse 
of notation we introduce variables S = ns/N, E = ue/N and / = nj/N which describe 
the population density of individuals in each disease state. Note that there is no need for 
a variable associated to the recovered state, since the conservation of total population 
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makes it a dependent variable. Applying the same system-size expansion as before (see 
Appendix A for details), we obtain the following effective SDE system: 

— = fi{l-S)-(3{t)SI + r]i{t), 
dE 

— = - (/i + a)E + r/2(t) , (38) 

^ = aE-(/i + 7)/ + r73(t), 

where 771.2,3 are Gaussian white noise variables with correlations 

1 

iV 



■ (39) 



/ S) + p{t)SI -iiE + (5{t)SI -fil \ 

B= -fiE + ^{t)SI ^{t)SI + {i2 + a)E -aE 

y —fil —aE aE + (/i + 7)/ j 

In this section we discuss the behaviour of the model in the deterministic limit 
—7- 00, before moving on to consider the full stochastic system in Section [321 When 
the infection rate is not seasonally forced, so P{t) = /3, there are a pair of fixed points. 
The first of these represents the extinction of the disease: S = 1,E = 0,1 = 0. The 
second fixed point has coordinates 

(a + /x)(7 + /x) fijl - S^) afijl-S^) 

^ = ' ^ = \ ' ^ = 7 — \ — u — < — ^ ' (^^'' 

ap a + jj, (a + /ij(7 + /i) 

and is referred to as the endemic state. We are concerned with the regime in which 
the endemic state is stable and the extinct state is unstable, which holds for a range of 
epidemiologically realistic parameter values. In general, the rate parameter /i controlling 
birth and death will be much smaller than the remaining rate parameters /3, a and 7, 
since this takes place on a much longer timescale than the disease dynamics. The 
parameter fi can therefore be utilized as an expansion parameter to simplify some of 
the expressions in the analysis. To first order in /i, the eigenvalues of the Jacobian at 
the endemic state are 

a{2a + (3) + 3a^ + 27^ 

Ai = - a + 7 -/i , 41 

[a + 7)"^ 

and the complex-conjugate pair 



0^^+^72 + 07(^+7) a(/3-7) .... 
^2,3 = /^ 7n — ^ — ^^ ±M//^ , ■ 42 

27(0 + 7)2 y a + 7 

Notice that we have a separation of time-scales: Re[Ai] ^ Re[A2,3], nieaning that Ai 
corresponds to a highly stable direction. This behaviour has been previously noted and 
exploited in the deterministic setting [22]. In addition, since /i is small, the imaginary 
parts of A2,3 are typically larger than the real parts, meaning that we may expect highly 
oscillatory trajectories in the neighbourhood of the endemic state. We can thus expect 
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0.0645 



Figure 3. Deterministic trajectories for the SEIR model. Left: No seasonal forcing 
{6 = 0), fixed point shown in red. Right: system in the presence of forcing, [S = 0.02), 
limit cycle highlighted in red. Remaining parameters in both plots are /3o — 1575 
year"-^, a = 35.84 year~^, 7 — 100 year"^ and fi = 0.02 year^^ which are standard 
choices for measles [T91I23]. 



the system to first collapse rapidly in the direction of the first eigenvector, followed by 
a slow, almost-planar, spiralling decay to the endemic state, as shown in figure |3l 

Introducing seasonal forcing of the infection rate creates an additional layer of 
complexity. A typical choice would be 

/3(t) = /3o(l + 5cos(27rt)) , (43) 

where Pq describes the basal infection rate, 6 is the forcing amplitude, and time is 
measured in years. The deterministic system will now not settle to the endemic state, 
but instead exhibit limit cycle behaviour. For the sake of simplicity, we will consider 
only parameter values which result in a single stable limit cycle of period T = 1 year. 

Similar to linear stability analysis of fixed points, there is a well-developed theory 
for analysing perturbations around limit cycles. Writing {S*{t), E*{t), I*{t)) for the limit 
cycle, we introduce the vector 

/ S{t)-S*{t) \ 



(44) 



E{t)-E*{t) 
\ /(t)-/*(t) J 

where the factor a/TV has been introduced for notational consistency with later sections. 
Since we are effectively linearizing the system, it should be noted that this factor has 
no effect on the following analysis. To first order then, the dynamics of ^ are governed 
by 

f = ^^^^^ ' ^^^^ 
where J is the time-dependent Jacobian of fl38l) . This equation may be solved using 
Floquet theory [24j. The key result of the theory states that solution trajectories 
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may be decomposed into the product of a periodic vector with an exponentially 
growing/ decaying amplitude. General solutions are then of the form 



(46) 



where n is the number of degrees of freedom in the model, q is a constant, qi{t) a 
periodic vector with the same period T as the limit cycle, and the value cTj determining 
the rate of growth/decay is referred to as the Floquet exponent. Akin to the eigenvalues 
of the fixed point Jacobian, Floquet exponents are indicative of the stability of the limit 
cycle in the time- varying directions qi{t)] perturbations to the trajectory will grow if 
Re[(Tj] > and decay if Re[crj] < 0. 

In all but the most trivial examples, obtaining the Floquet exponents and periodic 



vectors must be carried out numerically. The procedure (detailed in Appendix B) 
requires first computing a matrix whose columns are the independent solutions to ( l45l) . 
^{t) = (^i(^) ■ ■ '^nit)), and then determining the eigenvalues of XiT). If there is a 
rapid collapse along a stable direction, followed by slow dynamics in a subspace, then the 
columns of X{T) will be almost linearly dependent, since all solution trajectories quickly 
move towards the subspace. In turn, the real parts of eigenvalues of the matrix X{T) will 
differ by many orders of magnitude. Obtaining these eigenvalues accurately is crucial to 
the remaining analysis. However, if the disparity between the eigenvalues is too great, 
numerical procedures may not maintain sufficient accuracy in calculations involving the 
eigenvalues. This problem was previously highlighted in [19] within an analysis of the 
seasonally forced SEIR model. There, using the same epidemiological parameters as 
used in figure [3l the eigenvalues of X{T) differed by an order of 10^^. This prompted 
the authors to implement arbitrary precision numerical methods, at a considerable cost 
in computing time. We propose that the general method detailed in Section 12.31 can be 
used as a technique to remove this fast direction and hence circumnavigate the numerical 
difficulties encountered in the analysis. 



3.2. Stochastic treatment exploiting the slow manifold 

Beginning with the unforced case, we let Ai, A2 and A3 be as in fHT] H2|) and write 
Vi, V2, for the corresponding eigenvectors. We introduce the transformation matrix 
V = {vi + V3) i{v2 — I'a))^ and new variables 



(47) 



The Jacobian of the transformed system at the endemic fixed point takes the form 
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The nullcline for w is determined by manipulating fl29p into the form w = 6{z2, z^), with 
A copied from fl38|) . Although an explicit form for 6 can be found, the expression is too 
complicated to be worth reproducing here. 

To capture the effects of stochasticity, we introduce variables describing the 
fluctuations in the new coordinates, rescaled by a factor of \/iV, 

(49) 

Making this substitution in ( l38l) and keeping only first order terms in and fi, we 
find that ^ obeys 

§ = ^T+C(t), (50) 

where 

mXjit')) = Kt - t')B,, , z,j = 1, 2, 3. (51) 
The matrix B is given by 

B = VBV^ , (52) 

(S,E,I)=(S*,E*J*) 

where B is as in ( l39l) . Note that applying the constraint w = 9{z2,Z3) induces the 
relationship 

+ ^ = ^ (^^2 + 4^,^3 + 4^1 • (53) 



Expanding once more in large A^, we find 



f 06 - 39 

, dZ2 OZs 



(54) 

After elimination of the fast direction, ( 150|) becomes 

(55) 





where C2 and C3 now have correlation matrix B, which is related to B by 

We move on now to study the situation of seasonally forced infection rate. In 
principle, the calculations above apply only in the limit of small forcing amplitude 
(that is, 5 — )■ in (HH]) ). We learnt in the illustrative ecological example, however, 
that although our theory is developed to apply in the locality of a stable fixed point, 
it can continue to provide a useful approximation if this condition does not strictly 
hold. Applying that lesson to the present case, we modify ( 155|) to allow L and B to 
become functions of time, as dictated by the replacement (3 1— )■ P{t). Essentially, we are 
approximating the limit cycle by its projection on to the nullcline of the fast direction at 
the endemic fixed point of the unforced model, and then demanding that any stochastic 
fiuctuations remain on this nullcline. The application of such a coordinate change to 
the forced system can be further motivated by a consideration of the periodic matrix 
J{t) in (H5l) in the limit of small forcing. 
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A Floquet analysis of the deterministic part of the reduced system finds two complex 
conjugate Floquet multipliers, with the third disparate multiplier having been eliminated 
from the system. This system no longer suffers from the numerical difficulties which 
plague the full three-dimensional model. 

To quantify the effect of stochastic fluctuations in this model, we follow the standard 
procedure of computing the autocorrelation matrix C(r) of oscillations around the limit 
cycle, which has entries 

C{t)]^^ = (UtMt + r))dt- (56) 

Of course our reduced system ( l55l) is two-dimensional, meaning that the entries of C 
pertaining to must be deduced from (l54l) . The coordinate transformation applied 
at the start in P7|) may then be inverted to give the autocorrelation matrix for the 
fluctuations in S, E and /. 

The Fourier transform of the diagonal entries of the autocorrelation gives the power 
spectrum of oscillations, which provides a convenient visualization of the stochastic 
fluctuations. In figure H] we plot the power spectrum of fiuctuations in the number of 
infected individuals around the limit cycle, comparing between stochastic simulations 
and the theoretical prediction using the reduced-dimension model (|55l) . The agreement 
between the simulation and theory can be seen to be excellent; the spectra lie virtually 
coincident. The peaks in the theoretical spectrum are found at the same positions as 
those for the simulated spectrum. These are approximately given by 

., = 1±H!^, (57) 



where j is an integer and cxj are the Floquet exponents as defined in Appendix B 



This is in agreement with where peaks at these positions were also found. The 

overall benefit that was garnered by the procedure however was that of computational 
efficiency. The computation of the theoretical power spectra approximated from the 
reduced system takes only a small fraction of the computing time of the full system. 



4. Conclusion 



In this article we have introduced a systematic and general procedure to eliminate 
fast degrees of freedom in stochastic dynamical systems. The method was inspired 
by the highly successful theory of slow manifolds in deterministic systems, from which 
we borrow the basic notion of restricting attention to trajectories occupying a low- 
dimensional subspace. In the stochastic setting, we achieve this by enforcing the 
condition that the system state remains fixed to the nullcline of the fast direction. We 
apply this condition to SDE systems with correlated Gaussian white noise, obtaining 
a lower-dimensional system in which the fast degree of freedom has been eliminated. 
Importantly, the reduced system has the same basic form as the original, meaning that 
we have not inadvertently introduced any complicating factors such as non-Markovian 
processes or coloured noise. 
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frequency (year ) 

Figure 4. Power spectra for the fluctuations in the number of infected about 
the limit cycle. The analytical power spectrum calculated using the slow manifold 
approximation is plotted in red while the power spectrum from stochastic simulations 
is in blue. Agreement is such that the spectra are difficult to distinguish; the spectrum 
from simulated results is primarily discernible through its stochastic nature relative to 
the smooth analytical line. Dotted lines indicate the position of the peaks in the power 
spectra given by (f57| . Epidemiological parameters are /3o = 1575 year~^, a = 35.84 
year~^, 7 = 100 year^^ and /i = 0.02 year^^. The simulated spectrum is calculated 
as the average power spectrum of 1000 stochastic realizations, each lasting 200 years 
with a system size of iV = 10^. 

Although our procedure is by no means the only way to reduce the dimension 
of a stochastic dynamical system, it does offer some distinct advantages. We would 
argue that our approach is relatively simple: working in the SDE setting makes clear 
the analogy to deterministic slow manifold theory, and our basic idea (confining the 
stochastic system to the slow manifold) is intuitively easy to grasp. Moreover, the 
technique is generally applicable; it does not require an explicitly fast variable, since 
the fast direction is determined instead from a linear stability analysis. Nor does the 
technique itself require knowledge of the deterministic trajectories, as is necessary in [T5] . 
Finally, as we have shown, the approximation often remains successful even outside 
of the parameter regime in which it is developed. Confinement to the slow manifold 
does, however, remove any effects arising from the nature of the flow field away from 
the manifold. For example, in [15] it was shown that the way in which stochastically 
perturbed trajectories relax back can induce a small flow along the manifold; this effect is 
overlooked by our method. Looking to the future, there is the possibility to developing 
a general framework which includes effects of this type, and we hope that the work 
presented in this article will provide a significant step towards that goal. 
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Appendix A. 

In this appendix we give the details of the derivation of the SDE systems ([7]) and f l38|) . 
which are mesoscopic descriptions of the individual based models used as examples in 
sections 2 and 3, respectively. The general formulation of the method is described in |17j . 

Recall the model studied in section 2, with reactions listed in (|6]). The state of the 
system at a given time is specified by the vector n = {ux, ny), where nx and ny denote 
the number of individuals in populations X and Y, respectively. We write P{n, t) for 
the probability that the model is in state n at time t. To properly specify the model, 
we must write a master equation describing the time evolution of P. 

The reactions ([6]) define the possible ways in which the system can transition from 
the state n' to the state n; the rate with which these transitions occur is written T{n\n'). 
Specifically, we have 



The master equation is written in terms of the transition rates in the following way; 



While the master equation describes the dynamics of the system entirely, it is in general 
not analytically tractable. A particular realization of a process obeying the master 
equation may be simulated using the Gillespie algorithm |20]. While the simulation 
does not necessarily provide as deep an understanding of the system as an analytical 
treatment, it provides a useful comparison for the analytical techniques. 

To make analytical progress, we consider the limit of small e, applying a procedure 
known as the Kramers- Moyal expansion [TTlIlH]- Briefly, it involves the introduction of 
the rescaled state vector x = en, followed by a Taylor expansion of T and P around 
X. Truncating after two terms, results in a Fokker-Planck equation [26j describing the 
evolution of the probability density function; 



T{nx + l,ny|nx,ny) 
T{nx,nY + l\nx,nY) 
T{nx - l,ny|nx,ny) 
T{nx,nY - l|nx,ny) 



(1 - fl)nx + fiUy 

fiux + (1 - fJ')nY 

£(1/2 - p)nj^ + e(l/2 + p)nxnY 

e{l/2 + p)nxnY + e{l/2 - p)n^ 



(A.l) 




(A.2) 



n 



dP{x,t) 
di 



E d^. lM^)Pi^, t)] + -J2 [B.MP{^, t)] .(A.3) 
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The vector A{x) and the matrix B{x) can be calculated from the probability transition 
rates flA.l|) . It can be shown |T2] that the above equation is equivalent to the SDE 

^ = A{x) + ri{t), (A.4) 

defined in the sense of Itb [27j. Here r}{t) is as usual a Gaussian white noise term such 
that {r](t)) = and {'r]i{t)r]j{t')) = eBijS(t — t'). For the present example, the SDE 
formulation is given in ([7]). 

The procedure is much the same for the epidemiological model considered in section 
4. The transition rates are this time derived from the reactions fl37|) . and the state vector 
is now n = {ns,nE,nj). As is usual in this type of model, we choose (where N is 
the total number of individuals in the population) for the small parameter, e, used in 
the Kramers- Moyal expansion. (155]) is found as the result. 



Appendix B. 

The analogue of a linear stability analysis for systems with periodic components is 
known as Floquet theory [24] . It can also play an important role in the analysis of 
stochastic fiuctuations about a deterministic trajectory [19l[25]. In this appendix the 
general formulation of Floquet theory is discussed before the more detailed application 
to linear stochastic systems is given. 

Floquet theory gives the solutions to sets of linear differential equations in the form 
of (HSil . where J{t) is periodic with a period T. The general solution can be shown to 
be 

n 

^(t) = ^Qg,(t)e'^»*, (B.l) 

i=l 

where qi{t) is a periodic vector and cxj are termed the Floquet exponents of the system. 
The quantities pi = e°^*^ are termed the Floquet multipliers of the system. 

In particular one can work in a canonical form for calculational ease, with canonical 
quantities denoted with a superscript 0. The canonical form is constructed from n 
decomposed solutions to (l45l) such that ^f'\t) = qf'\t)e'^'^. A fundamental matrix of 
these solutions may then be introduced along with matrices Y^'^^ and Q^^\ For the case 



n = 3 these may be expressed as 

X(°) = [|f)(t),^f(t),^f(t)], (B.2) 

Xio) = Q(o)y{o)^ (B.3) 

Q(°) = [gr(t),qf(t),gf(t)], (B.4) 

= Diag[e^^% (B.5) 



A method for obtaining the Floquet multipliers /ij along with the canonical form of 
the solutions is now required. Obtaining both is dependent on the determination of a 
matrix known as the monodromy matrix, which we shall now discuss. 
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The monodromy matrix, D, is defined such that X{t + T) = X{t)D, for any 
fundamental matrix X{t) constructed from hnearly independent solutions to fj45|) . It can 
be shown that while the monodromy matrix is dependent on the fundamental matrix 
chosen, its eigenvalues are not [23]. The eigenvalues of D are p,, the Floquet multipliers 
of the system. Further, if a matrix W is constructed from the eigenvectors of D, the 
canonical fundamental matrix X^^'^ (t) is related to a general fundamental matrix via 
X^^\t) = X{t)W. Therefore, the monodromy matrix allows the canonical fundamental 
matrix X*^°-'(t) to be determined from a general fundamental matrix X(t), along with 
the matrix Y^^\ From these the periodic matrix Q^'^\t) can also be deduced. 

If the system under consideration displays limit cycle behaviour, the preceding 
analysis must evidently be carried out numerically. In general a fundamental matrix 
obtained numerically will have to be transformed into canonical form by a numerical 
determination of the monodromy matrix, D = X~^{t)X{t + T). For a system with 
initial conditions t = 0, X{0) = I, this simplifies to = X{T). 

Now the stochastic system can be considered; 

^ = J(t)| + ry(t) , (B.6) 



where r){t) is a vector of Gaussian white noise terms defined as in Appendix A[ except 



that now the noise covariance matrix depends explicitly on time through the varying 
parameter /3(t); {i]i(t)rij(t')) = eBij{t)5{t — t'). The solution may be constructed as a 
sum of the general solution to f l45|) along with a particular solution, so that 

t 

^(t)=X(°)(t)^W+X(°)(t) j [X^''\s)y\{s)ds, (B.7) 

to 

or, setting the initial conditions in the infinite past and making a change of integration 
variable s ^ s' = t — s 

t 

i{t) = Q^^\t) j F(°)(s') [Q^^'Kt - s')Y\{t - s')ds' . (B.8) 

to 

In the course of the analysis conducted in Section [31 ^{t) represents some stochastic 
fluctuation around limit cycle behaviour. An obvious quantity of relevance is the power 
spectrum of such fluctuations. To obtain the power spectrum, one first calculates the 
two-time correlation function C(t+r, t) = (^(t+T)^"^(t)); substituting (]B.8|) one obtains 

1 T 



with 



a,{t + T,t) = g(°)(t + r)r(°)(T)A(t) [Q^'\t)] , (B.9) 



A(t)= Y^'^\s)T{t-s)Y^'^\s)ds (B.IO) 



to 



and 



T{s) = [Q^'\s)]-' Bis)\[Q^'\s)YX ■ (B.ll) 
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The correlation function, C(r) is then simply related to the two-time correlation function 
by 

C(r) = l J Cit + T,t)dt. (B.12) 



In turn, the Weiner-Khinchin theorem tells us that the power spectrum, P{uj), is simply 
the Fourier transform of the correlation function, and so 

Pi{uj) = j Cu{r)e'''^dT. (B.13) 

The intermediate steps are left to the reader, but full details are found in [28j. A key 
point to note is that Equations (IB.THB.lT]) hold only for the canonical matrices X^^\ 
Q(o) and 
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